Quasiparticle energies for large molecules: a tight-binding GW approach 
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We present a tight-binding based GW approach for the calculation of quasiparticle energy levels 
in confined systems such as molecules. Key quantities in the GW formalism like the microscopic 
dielectric function or the screened Coulomb interaction are expressed in a minimal basis of spherically 
averaged atomic orbitals. All necessary integrals are either precalculated or approximated without 
resorting to empirical data. The method is validated against first principles results for benzene and 
anthracene, where good agreement is found for levels close to the frontier orbitals. Further, the size 
dependence of the quasiparticle gap is studied for conformers of the polyacenes (C4n+2-ff2n+4) up 
to n = 30. 



I. INTRODUCTION 

Density functional theory (DFT) has nowadays be- 
come the standard tool for the description of ground state 
properties of such different systems as atoms, molecules, 
clusters or bulk materials. Part of the success stems from 
the fact that in the DFT electron correlation is in princi- 
ple exactly covered at the level of an effective one-particle 
Hamiltonian. Difficulties arise however, when the orbital 
energies obtained from the solution of the Kohn-Sham 
(KS) eigenvalue problem are interpreted as approximate 
quasiparticle energies, i.e., the energies associated with 
addition or removal of an electron. A well known ex- 
ample is the severe underestimation of the band gap in 
insulating solids or semiconductor crystals . The same 
problem is also present in molecules. Although it has 
been shown that Koopmans theorem also holds for the 
highest occupied molecular orbital (HOMO) in DFT p|, 
ionization potentials come out too small in practical cal- 
culations. This fact has been traced back to the wrong 
asymptotic behavior of common approximate exchange- 
correlation (XC) functional like the local density approx- 
imation (LE)A) However, the KS gap can be shown 
to represent a first-order approximation to optical exci- 
tation energies 3, and thus the KS gap will be different 
from the quasiparticle gap even if the exact XC functional 
is used. 

As an alternative to DFT, many body perturbation 
theory in the approximation of Hedin j5| has been ex- 
tremely successful in the prediction of quasiparticle spec- 
tra. In this so-called GW method the description of 
the electron-electron interaction is approached in a dif- 
ferent way than in DFT. While the exchange energy is 
exactly calculated like in Hartree-Fock theory, correlation 
is accounted for by an energy dependent dielectric func- 
tion which reduces the coulomb interaction between elec- 
trons. The scheme is nearly self- interaction free and pro- 
vides asymptotically correct potentials. Consequently, 
the band gap problem of the DFT is absent in the GW 
and also ionization potentials and electron affinities of 



molecules are computed in good accord with experiment 



In conjunction with the solution of the Bethe-Salpeter 
(BS) equation the GW approximation can 

also be used to calculate charge neutral quasiparticle ex- 
citations , i.e., optical spectra. In this context, the time- 
dependent generalization of DFT ,1^] , which in principle 
provides the same information as GW-I-BS has in practice 
been found to give inferior results for the absorption spec- 
trum of solids 14j. Also for molecules, time dependent 
DFT fails in the description of charge transfer excitations 
[m , which was attributed to the locality of common XC- 
functionals [IE 113. El- Again, GW+BS which contains 
the correct long range electron-hole interaction should be 
able to remedy the problem. 

Another example where DFT orbital energies are 
widely used but lead to problematic results is given by 
transport calculations in molecular devices. Here the cal- 
culated currents differ typically by orders of magnitudes 
from the experimentally found values [T^ l20l |21| . Since 
the current-voltage characteristics of such systems de- 
pend critically on the HOMO-LUMO gap, a major im- 
provement is expected when GW quasiparticle energies 
would be used instead of DFT orbital energies. 

It is clear from the forgoing discussion, that although 
the GW approximation was developed in the context 
of solid state theory, its application to molecular sys- 
tems is becoming more and more important. In fact, 
implementations for systems with translational symme- 
try using plane wave basis sets can directly be applied 
to finite systems when very large super-cells are used. 
Nonetheless, the use of localized atomic-like basis sets 
is clearly more adapted to the problem and such imple- 
mentations have been utilized quite successful in the past 
years 6. 7. 22.. .23j. But even with this improvement, the 
numerical complexity of the GW equations limits a first 
principles evaluation to rather small system sizes of tens 
of atoms. So applications like transport calculations for 
molecules, where a sizable amount of the atoms com- 
prising the leads need to be included, or molecular dy- 
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namics in the excited state are currently not feasible. It 
would therefore be desirable to have an approximate GW 
scheme which nevertheless captures the essential physics 
of the underlying theory. The purpose of this paper is to 
propose such a scheme. 

In the past years a number of different approximated 
GW methods have been introduced and successfully ap- 
plied 5, 24, 25, 26, 27, 28, 29, 30,, 31, 32]. The main 
difference to this earlier work is that we focus on molec- 
ular applications with a real space implementation. Fur- 
thermore we try to avoid empirical parameters in order 
to achieve a higher transferability. 

Like GW calculations usually employ DFT energies 
and wavefunctions as zero order approximations to the 
quasiparticle quantities, our approach is based on the 
Den sity Functional based Tight- Binding (DFTB) method 
[3^ l34|. which itself is an approximation to DFT. The 
DFTB scheme has been shown to provide reliable results 
on a variety of systems classes ranging from molecules 
to solids at a highly reduced numerical cost compared 
to DFT calculations. In section^ the DFTB method is 
briefly introduced together with a detailed description of 
the approximations in the various quantities involved in 
the GW formalism. The accuracy and main shortcom- 
ings of our approach are then examined in section IIIII 
where it is applied to a prototype series of 7r-bonding 
molecules, the polyacenes. 



II. METHODOLOGY 

The GW method has been extensively discussed in the 
literature and several reviews are available 14, 36, 37, 
Is^ls^- The main goal is to solve the Dyson equation: 

{Ho + ner))\^^?') = erH?% (1) 

for the quasiparticle energies ef^ and wavefunctions 

Here Hq is the Hartree Hamiltonian and the so- 
called self-energy E is a nonlocal and energy dependent 
operator, which accounts for exchange and correlation 
effects. It thus can be seen as a replacement for the lo- 
cal exchange-correlation potential v^c in the DFT Kohn- 
Sham (KS) equations. In the GW approximation of 
Hedin, the self-energy is given as a product of the single- 
particle Greens function G and the screened Coulomb in- 
teraction W. As these quantities, as well as the Hartree 
Hamiltonian, depend on the quasiparticle wavefunctions, 
the Dyson equation [Eq. (^] has to be iterated until self- 
consistency is achieved. However, since the DFT one- 
particle wavefunctions are usually very similar to the fi- 
nal quasiparticle ones and moreover self-consistency not 
necessarily improves the result [ioj . Eq. may be sim- 
plified to 
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where V'i are KS orbitals, non-diagonal elements of the 
Dyson Hamiltonian in the basis of DFT states are ne- 



glected and the energy dependence of the self-energy has 
partially been accounted for by the renormalization con- 
stant Z,: 



duj 



(3) 



In our approach Eq. (0) is now subjected to further 
approximations which are presented separately for each 
term in the following. 



A. The DFT orbital energies ef 

As already mentioned, the Kohn-Sham energies ef^^ 
and orbitals are obtained from the DFTB method, which 
has been presented in detail earlier (for a review see 
Ref. ^3). Here we only describe the method to the ex- 
tent necessary to motivate the remaining approximations 
in this work. In the DFTB, the Kohn-Sham states tpf^^^ 
are expanded in a linear combination of atom centered 
orbitals 6,, : 



(4) 



which are obtained from a preceding DFT calculation on 
neutral atoms. Here jj. := {Aim} is a compound index in- 
dicating the atom on which the basisfunction is centered, 
its angular momentum I and magnetic quantum number 
TO. Later, also quantities which depend only on A and 
I appear. These will be denoted with a corresponding 
index p, := {Al} throughout the paper. 

Since atomic orbitals are usually too long ranged to be 
used directly in a molecular calculation, the atomic DFT 
Hamiltonian is augmented with a confining square po- 
tential to compress the wavefunction outside of a given 
radius tq (usually twice the covalent radius of the ele- 
ment), while ensuring the desired cusp conditions inside 
[3^ . From the atomic valence states a minimal basis of 
s and p orbitals is then chosen, although also d-orbitals 
are included when necessary, e.g. for second row elements 
. With the help of the expansion the Kohn-Sham 
equations of DFT can be written: 



sec 



(5a) 

(5b) 
(5c) 



where the overlap matrix S^i, has been introduced and 
the Hamiltonian is divided in two parts. The first part 



H^j^ is approximated as follows: 



^frcc atom 

= \ {M^)\HDFT[p'A + P%]\Mr)) 



AieA,j/eB (6) 
otherwise 
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The KS Hamiltonian in Eq. © contains as usual the 
kinetic energy, the electron-nuclei attraction and the 
Hartree as well as exchange-correlation potential and de- 
pends only on the atomic densities of atoms A and B. 
This means, that besides the crystal field terms also all 
three-center terms are neglected. The onsite elements are 
given as atomic orbital energies obtained from a DFT cal- 
culation without confining potential to ensure the right 
dissociation limit. The integrals in Eq. © are numer- 
ically evaluated and tabulated for varying distance be- 
tween atoms A and B. 

The second part of the Hamiltonian ifKcjl corrects for 
the fact that the molecular density differs from a simple 
superposition of atomic densities. In order to estimate 
this difference, spherical averages over basis functions be- 
longing to one angular momentum shell are build 



FAiir) = 
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and used in a Mulliken type approximation 

</>^(r)0,(r) « i V {F-,{r) + F,{r)) , 



(7) 



(8) 



to represent the molecular density. The latter is then 
constructed from point charges 



i ft 
I 

with qp = ^ y^^c^iC^.iSf_,^; 



(9) 
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an expansion, which despite of its simplicity takes the 
different spatial localization of e.g s and p-orbitals into 
account. Based on these considerations, the difference 
between the true molecular density and superimposed 
atomic densities can be estimated with net Mulliken 
charges Aq^ = qj, 
term .35.1 : 



^atom g^j^j leads to the correction 



ttSCC 



1 



'S'M^y^(7u5 + 7pj)Agj, 



S 



(10) 



as shown in more detail in Ref. p34l |. 

The term 7 describes the interaction of two electrons 
in the orbitals ft on atom A and v on atom B, including 
the effects of exchange and correlation: 



dp{r') 



Jf,,{\JlA~RB\,U"U^), 



Fp(r')drdr' 

(11a) 
(lib) 



and is approximated by considering two limiting cases. 
For large distances between the two atoms, the integral 
(|lla|l simplifies to a pure Coulomb interaction of two 



point charges, since the v^c contribution dies off rapidly. 
For short distance on the other hand, Eq. Hlla|) becomes 
an atomic integral Uh, which can be easily calculated 
numerically for each element. From these limiting cases 
a simple interpolation formula was derived in Ref. |34j . 
which is a function of the atomic parameters and 
the interatomic distance only. Since the Mulliken net 
charges Aqg depend on the molecular orbital coefficients 
Cf^i, Eq. H5a|l has to be iterated until self-consistency. As 
a result the orbital energies ef^"^^ needed in Eq. (O are 
obtained. 



B. The self-energy Si(e) 



We calculate the self-energy in the GW approximation 



by: 

S(r,r',e) = ^ J e-°+Go(r, r', e - c.) Wir,r',u;)dco, 

(12) 

where Go is the single particle Greens function built from 
DFTB wavefunctions and W — e~^v is the screened 
Coulomb interaction, while v is the bare one. For the 
following it is beneficial to divide the self-energy into two 
parts as S = iGov + iGo{e~^ — l)v. The latter term 
denoted E'^ is energy dependent and describes dynamical 
correlation effects, while the former term provides the 
major part of the self-energy. For E^ the frequency inte- 
gration in Eq. (|12|1 can be carried out easily and yields 
in the KS basis the usual Hartree- Fock exchange energy 
for orbital i: 



E 



drdr' . 



(13) 



Since in contrast to empirical tight-binding schemes the 
basis functions (f>f^ are available in the DFTB method, 
one could in principle calculate Eq. (|13|l directly from the 
wavefunctions. In this way the method would scale like 
first principles schemes with N'^, where N is the number 
of basis functions. Therefore we seek for an approximate 
solution and note that after expansion of the KS states in 
atomic orbitals, Eq. (|13|l contains products of basis func- 
tions which are in general located on different atomic cen- 
ters. An important simplification can thus be achieved, 
when the Mulliken approximation [Eq. (jHJ] is applied to 
the integral. Introducing the following notation for the 
matrix elements of a general two-point function in the 
basis of squared and spherically averaged DFTB atomic 
orbitals: 

[fU = // F^(r)/(r,r')F,(r')drdr', (14) 



we then arrive at the following simplified expression for 
Ef: 



(15) 
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Here the q^^ are generalized MuUiken charges 
1 ' 

which provide a point charge representation of the over- 
lap between two molecular orbitals i and j. The impor- 
tant observation is now that the matrix of the Coulomb 
interaction is equal to the definition of the 7-functional 
in Eq. (|lla|l . when the contributions stemming from the 
XC functional are removed. In other words, the func- 
tional form of 7 can also be used for [u] , if the atomic pa- 
rameters are replaced by the parameters IJ^^ which 
incorporate only the classical Coulomb interaction. We 
calculate these electron repulsion integrals directly from 
the DFTB basis functions using the algorithms presented 
in Ref . 43] . The parameters for each angular momentum 
are set to an average over the integrals for different com- 
binations of the magnetic quantum numbers belonging 
to that shell. 

The main drawback of the MuUiken approximation in 
Eq. © is, that onsite integrals of the exchange type are 
completely neglected. These, however, contribute around 
10 % to the final exchange energy. Similar to the pro- 
ceeding in the quantum chemical INDO approach [4J| we 
therefore include all non- vanishing onsite integrals, lead- 
ing to the following final form for Ef |45l |: 



3 A"^ 



1.1 iLu 



(17) 



While in the INDO approach the necessary integrals are 
taken as empirical fitting parameters, we compute them 
from the atomic basis functions. More precisely, the pa- 
rameters are calculated from an uncompressed wavefunc- 
tion in order to be consistent with the onsite definition 
of the DFTB Hamiltonian matrix elements. The values 
used in this study are given in Tab. U together with the 
and parameters. 
Let us now turn to the correlation contribution of the 
self-energy E^, which is much harder to evaluate, since 
it amounts to a multi step procedure. First we construct 
matrix elements of the electronic polarizability in the 
random-phase approximation according to: 

occ virt / \ / 

k I \ a / \ B 



(18) 



^DFTB ,DFTB 
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where we again used the MuUiken approximation of 
Eq. ijHJl and introduced the overlap matrix of spherically 



TABLE I: The atomic electron-electron interaction inte- 
grals Ul^ and f/f^, as well as the exchange integrals 
{4>im4'i'm'\4>im<t>i'm') used in this study. Results are given for 
free and compressed atomic basisfunctions, as defined by the 
confinement radius ro (see Sec. Ill Al and Ref. [s^]). The same 
compression is used in the calculation of the Hamiltonian and 
overlap matrix elements. 



Element 



Parameter 



ro [a.u. 



Value [eV] 



Hydrogen 



Carbon 



ur 

(t>004>lm'\4^004>lm') 
4>lm' \4>\m4>l m 



00 
00 

3.0 

00 

CX3 
00 
00 

2.7 
2.7 

00 
00 



11.06 
15.39 
21.36 

10.81 
10.81 
15.66 
14.15 
17.98 
18.72 
3.01 
0.75 



averaged DFTB basis functions S'^p — J Ffi{r)Fp{r)dr, 
not to be confused with the overlap appearing in the KS 
equations (|5b|l . 

The quantity S never needs to be constructed, since it 
appears only in intermediate quantities and falls out in 
the final equation for the screened Coulomb interaction 
we are aiming at. 

In a next step we obtain the dielectric function in ma- 
trix notation as: 



(19) 



For systems with translational symmetry the dielectric 
matrix is hcrmitian in reciprocal space. This fact is used 
in plasmon-polc models 24, 46, 47J to simplify the fre- 
quency integration in Eq. (|12|l . which is numerically de- 
manding due to the complicated pole structure of 
along the real axis. In these models the inverse of the 
dielectric matrix is diagonalized and the eigenvalues are 
assumed to be a simple function of the frequency, while 
the eigenvectors are frequency independent. Free param- 
eters of the model are either obtained from sum rules or 
by diagonalizing at different test frequencies. It is 
then easy to perform the frequency integration analyti- 
cally to obtain the self-energy. 

However, in the present real-space approach the in- 
verse dielectric matrix is not symmetric. In Ref. [23| this 
problem was circumvented by introducing an auxiliary 
symmetrized dielectric matrix, while we proceed by not- 
ing that the screened Coulomb interaction W: 

[Wiu;)]^SHc.)]-'[vl (20) 

has the desired property of being symmetric. Applying 
the plasmon-pole approximation to [W — v], we finally 
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arrive at the following expression for the correlation con- 
tribution to the self-energy for orbital i: 



n s \ fj. I 



1 



1 

,DFTtt_ 



n € occ 



n G virt, 



(21) 



where $ denotes the eigenvectors of [W — v] , while zj 
and Ug are the mentioned parameters of the plasmon- 
pole model, as defined in Ref. They are determined 
by diagonalization of [W — v] at zero frequency and one 
frequency on the imaginary axis. We checked that the 
actual values chosen have little impact (< 0.1 eV) on the 
final quasiparticle energies. 

Based on the self-energy, the renormalization constant 
Zi from Eq. |2Jl is then obtained by a simple numerical 
differentiation. For the molecular structures we studied, 
Zi is usually roughly 0.85, which is close to the values 
reported for bulk systems |2J|. However, for certain un- 
bound virtual orbitals, Zi can decrease to as much as 
0.5. 



C. The exchange-correlation contribution vf^ 

We complete the description of our method with an 
investigation of the contributions to the quasiparticle en- 
ergies arising from the exchange-correlation potential, de- 
noted vf^lpv]. As indicated, w^c is evaluated at the va- 
lence density pv , consistent with the fact that the summa- 
tion in the exchange energy [Eq. H15|l ] is carried out over 
valence orbitals j only. As pointed out in Ref. 0i the 
core contribution to the exchange energy is not neglible 
and this holds also for the core contribution of w^'^. How- 
ever, even for exchange-correlation potentials commonly 
used today, which are far from exact, both core contri- 
butions cancel to a large degree when computing quasi- 
particle energies according to Eq. 

In analogy to the derivation of the DFTB method, we 
now expand w^*^ around the density pJJ , which is a super- 
position of atomic valence densities. With p^, = p^, + Sp 
we obtain: 



vriPv]^ J mr)\' vr[p°{r)] dr + 
2 Svrip^ir)] 



5pv{r 

^c^,^c^,vll[pl]+"^q^ 



j^6p{r') drdr' + 0{6p^) (22a) 



}1V 



Sp 



Aqi, (22b) 



In going from Eq. (|22af) to Eq. 122bl) . the Mulhken ap- 
proximation was again employed and matrix elements of 
the exchange-correlation kernel Sv^'^/Sp were introduced 
in the notation of Eq. (djl. The first term of Eq. Ij22b|) 
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FIG. 1: The electron-electron interaction integrals for hydro- 
gen as a function of distance. Shown are results including 
both Coulomb and exchange-correlation interactions (7^), 
pure Coulomb (7'^'^) and the negative of the pure exchange- 



correlation interaction ( 
Eq. CT . 



-Sy'^'^/Sp = 7*"^ — 7 ): according to 



is now exactly treated like the Hamiltonian in the DFTB 
scheme. That is, only the two-center terms are kept and 
the onsite values are calculated from uncompressed basis 
functions and atomic densities. Then the integrals are 
numerically evaluated and tabulated in the usual Slater- 
Koster form as a function of interatomic distance. 

The second term in Eq. Ij22b|) is the counterpart of 
H^^'^ in Eq. ^U^. If we set 



~5i 



= 7^p(|RA-Ri3|,C/f,C/if) 



Rf 



(23) 



the long range 1/R tail of the two 7-functions cancels 
(see Fig. ^ , and one is left with a short ranged represen- 
tation of the exchange-correlation kernel without intro- 
ducing any new parameters or integral approximations. 
Moreover the w^'^ contribution of the self-energy now can- 
cels all related contributions in the orbital energies ef^"^^ 
as it should be 49]. 

At this point, all the necessary ingredients to cal- 
culate quasiparticle energies within the DFTB scheme 
have been presented. In the next section we analyze the 
strengths and weaknesses of the method by applying it 
to the polyacene series. 



III. APPLICATIONS 
A. Comparison to first principles results 

The polyacenes (C4„+2-ff2n+4) are linear chains of anel- 
lated polycyclic aromatic hydrocarbons, as shown in Fig. 
121 The monomer (n = 1) is benzene. Naphtalene is with 
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TABLE II: The difFerent contributions to the quasiparticle (QP) energies and the QP energies themselves obtained from the 
method described in this work (DFTB) as well as from first principles calculations using Gaussian type orbitals (GTO). Shown 
are results for some levels close to the frontier orbitals of the benzene and anthracene molecules. All energies in eV. 
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FIG. 2: Schematic viewgraph of the polyacenes: n is the num- 
ber of monomers. 



n — 2, anthracene n = 3, tetracene n = 4 and so on. 
These systems have received much attention because of 
their potential use in efficient organic thin film devices. 
Theoretically t hey have been characterized quite widely 
[sol ISll Isa . 153 . Isl . Issl ISq and recently also an investiga- 
tion in the context of the GW approximation appeared, 
where polymorphs of the pentacene crystal were analyzed 
in terms of their optical spectra [s^ . 

Here, the polyacenes are chosen as a prototypical tt- 
system to explore the accuracy of our approximations. 
To this end we fully optimized the different structures 
at the DFTB level of theory without imposing symme- 
try constraints. The obtained geometries are in excellent 
agreement with a recent DFT study on the polyacenes 
from n = 1 to 5, in which the hybrid functional B3LYP 
and the accurate 6-31 IG** basis set was employed [s^ . 
For all the molecules studied, we find a mean deviation 
in the bond lengths of no more than 0.005 A. Then the 
quasiparticle spectrum was obtained according to the ap- 
proximations in the last section. For comparison, also 
first principles GW calculations in the Gaussian type or- 
bital (GTO) implementation of Ref. were carried out. 
In the ab-initio GW calculation, the wave functions have 



been represented by s and p Gaussian orbitals on carbon 
(decay constants: 0.12, 0.4, 1.0, and 2.8 atomic units, 
i.e., 16 orbitals per atom) and by s and p orbitals on 
hydrogen (decay constants 0.15, 0.4 and 1.0, i.e., 12 or- 
bitals per atom). The two-point functions occurring in 
the GW scheme are represented by s, p, d, and s* or- 
bitals on carbon (decay constants 0.2, 0.6, 1.6, and 4.0, 
i.e., 40 orbitals per atom) and on hydrogen, as well (de- 
cay constants 0.25 and 0.7, i.e., 20 orbitals per atom). In 
both the DFTB and first principles calculations the LDA 
exchange-correlation functional was used. 

The results for benzene and anthracene are listed in 
Tab. ^for a small number of states around the frontier 
orbitals, according to the energy partioning of Eq. jSjl. 
Focusing first on the orbital energies, we find a very good 
agreement between the DFTB and the first principles re- 
sults. This might be surprising considering the limited 
basis set employed in the former approach. However, 
the DFTB basis consists of optimized atomic orbitals 
rather than simple Gaussian type orbitals. Moreover, 
the approximations underlying the DFTB method seem 
to be justified due to a stable error cancellation for the tt- 
orbitals. This hold true to a lesser extent for cr-orbitals, 
like the i?2g state in benzene, where an error up to 0.6 eV 
is found. Not unexpected, difficulties are also observed 
for unbound virtual orbitals like the B2g state in benzene, 
since the description of the continuum is very sensitive 
to the quality of a finite basis set. 

Next, we turn to the exchange-correlation energy per 
orbital. Here we find in comparison with the first princi- 
ples results, that the DFTB values are in general too pos- 
itive by roughly 10 %. However, the error even reaches 
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20 % for the cr-orbitals of benzene. We attribute this fail- 
ure to the neglect of crystal field terms in our approach, 
which is likely to have different effects on orbitals of tr 
and TT symmetry. In fact, we calculated elements of the 
type {4>^\ Vxc{pb) \'t'v) ^-i^d found, that integrals where 
A represents a hydrogen atom and B a carbon atom are 
significantly larger than in the reversed situation or inte- 
grals where both A and B stand for carbon atoms. As the 
latter two types of integrals occur in the calculation of 
TT-orbitals of the polyacenes, while the first one is impor- 
tant for cr-orbitals, the missing of the crystal field terms 
is likely to be the source of error here. 

Considering now the exchange contribution to the self- 
energy S^, similar trends are found. Compared to the 
ab initio results, the DFTB values are slightly too pos- 
itive. Since the terms and w*^ contribute in Eq. Q 
with opposite signs to the quasiparticle energies, a sta- 
ble error cancellation is expected. A larger deviation is 
found again for the E2g state in benzene, where an error 
up to 4 eV occurs. Since the exchange integrals depend 
strongly on the atomic repulsion integrals Uee in our ap- 
proximation, the error could be reduced by enlarging this 
parameter for hydrogen without loosing the good perfor- 
mance for the TT-orbitals. However, we hesitate to treat 
the Uee values as empirical parameters, because of loss of 
transferability. Instead, one should look for a better ap- 
proximation of the two-electron integrals. In our approx- 
imation the density 4>^{r)(p,y{r) is represented by a su- 
perposition of spherical charge densities. Consequently, 
the two-electron integrals are given by simple monopole- 
monopole interactions, thus neglecting any angular mo- 
mentum dependence. A natural next step would be to 
include higher order terms in a multipole expansion of 
the density <j>f^{r)(p^{r), as it is done in the semiempirical 
MNDO method developed by Dewar and Thiel [5a |. 

Next, the final quasiparticle energies are discussed. For 
the TT-orbitals the mean deviation of the DFTB results 
from the ab initio values is 0.4 eV, with errors decreasing 
when the system size is increased. As could be already 
expected from the forgoing discussion, the description of 
cr-orbitals is less satisfactory in the current state of ap- 
proximations. For the E2g orbital of benzene an error of 
1.25 eV is observed. For the unoccupied levels however, a 
very nice agreement between first principles and DFTB 
results is evident. Clearly, this is a consequence of an 
error cancellation between all terms in Eq. ||2Jl, since e.g. 
the correlation contribution is systematically under- 
estimated in the DFTB scheme. 

In this context it is interesting to investigate if a more 
advanced treatment of the Dyson equation Q leads to 
better results. In fact, it has been found that the asso- 
ciated wavefunctions of orbitals which are bound at the 
DFT level of theory, but unbound at the QP level, dif- 
fer considerably. This is in contrast to the assumptions 
made in the derivation of Eq. (0) from Eq. and hence 
the full QP Hamiltonian needs to be diagonalized in these 
cases and self-consistency with respect to the energy de- 
pendence of T, must be achieved. Following this approach 
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FIG. 3: Quasiparticle (•) and DFT gap (a) for the lowest 
energy conformer of the polyacenes as obtained in the DFTB 
approximation. Lines are guides to the eyes. Shown is also 
the difference of experimentally determined electron affinities 
and vertical ionization potentials (o) from Ref. |60ll . where 
they were available. 



earlier investigations of this point reported shifts of the 
LUMO level up to 0.8 eV [ll Hi]. We also performed 
such calculations for benzene and found that even for 
the E2u state the diagonalization changes the QP spec- 
trum by less than 0.01 eV. This can be understood as as 
consequence of the minimal basis set we employ, which 
does not provide enough flexibility to describe the relax- 
ation towards delocalized states. Considering the energy 
dependence of the self-energy, it can be stated that the 
approximate treatment of Eq. ^ using the renormaliza- 
tion constant Z is quite successful, as we find deviations 
less than 0.2 eV from the self-consistent solution of the 
Dyson equation. 



B. Size dependence of the quasiparticle gap 

After validation of the method, wc now turn to a first 
application and analyze in the following the quasiparticle 
gap = Clumo ~ ^homo ^ function of chain length. 
The first observation which can be drawn from Fig.|3is 
that the DFTB quasiparticle gap is in very nice agree- 
ment with the experimental data, which provides some 
confidence that the general trends we are looking for are 
correctly described. Furthermore, Fig. 13 shows that the 
DFT gap is continously decreasing and almost vanishes 
for n = 19 monomers. As the length increases, the geom- 
etry of the innermost part of the chain resembles more 
and more that of two coupled polyacetylene chains with 
equal bond lengths, as schematically depicted in Fig. 0] 
The vanishing of the DFT gap can thus be understood 
in terms of a simple particle-in-a-box model. 

In stark contrast to the DFT gap, e^^p remains finite 
for increasing chain length and it seems worthwhile to 
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a) 
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FIG. 4: Schematic representation of the lowest energy confor- 
mations of the polyacenes found in this study. The aromatic 
structure a) is the most stable for n < 19, while the Peierls 
(Z)-distorted structure b) is energetically favoured for n > 20. 



explore the physical origin of this different behaviour in 
a short digression. In fact, QP energies can be directly 
compared to results from photoemission and inverse pho- 
toemission measurements, i.e., they include the effect of 
an extra particle, while DFT is a pure N-electron theory. 
Delerue et al. pointed out that for nanocrystals the size 
dependence of the difference A = e^^^p — e^p-^ can be es- 
timated on the basis of classical electrostatic arguments 
|3l|. Considering the interaction between the extra par- 
ticle and its induced surface charge on the nanocrystal, 
they arrived at the following formula: 



1 



1 



e2 e2 fe(R)-l 

R e{R)R \£iR) + l 



(24) 

where e{R) is an effective dielectric constant and A;, is 
the bulk value of A. In order to apply Eq. H18() to the 
polyacenes, we took R to be half of the chain length 
and obtained e(i?) by averaging the microscopic dielec- 
tric function in Eq. I|19|) . The obtained values increase 
from 1.72 for n = 1 to 2.14 for n — 19, which reflects the 
decreasing band gap. A fit of Eq. ^| to our QP results 
is shown in Fig. [S] and leads to a value of 2.18 eV for 
A;,. Taking into account that the DFT gap is vanishing 
for n —> oo, we therefore predict a QP gap of the same 
value for an infinite chain in the aromatic structure of 
Fig. 21 Inspection of Fig. O further reveals that for n > 
4 the agreement between Delerue 's formula and the QP 
results is excellent. The fact that Eq. H18|) . which was 
developed in the context of nanocrystals also holds for a 
quasi one-dimensional system like the polyacenes is quite 
remarkable. 

We now continue the discussion of Fig. 13 Between 
n = 19 and n = 20 HOMO and LUMO cross, which 
has important implications for the geometrical as well as 
electronic structure of the system. Remaining in the pic- 
ture of polyacene as coupled polyacetylene, we observe 
that the equal C-C bond lengths found for n < 20 turn 
into alternating single and double bonds for larger n as 
depicted in Fig. ^ i.e., the system undergoes a Peierls 
distortion. In contrast to polyacetylene, where the bond 
alternation is found to be around 0.08 A |£1,], the effect is 
much weaker here with a value of less than 0.008 A. Nev- 
ertheless, the dimerization leads like in polyacetylene to 
an opening of the DFT gap, which tends towards a small 
but finite value for the infinite chain. Inspection of Fig.|2| 
shows that also the quasiparticle gap differs significantly 
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FIG. 5: Fit of Eq. JTHl (x) to the difference between QP and 
DFT gap as obtained in this work (•). 



for the distorted and undistorted structure. This fact 
could be used in experiment to discriminate between both 
polymorphs, since we find in line with the MP 2 results 
of Cioslowski '52'|, that the two forms are energetically 
quite close and should therefore coexist in real samples. 
We should mention however, that up to now only poly- 
acenes up to n = 6 could be isolated, since larger chains 
are highly vulnerable to photooxidation ■ 

The results of this section can be summarized as fol- 
lows. First, the aromatic from of the infinite polyacene is 
predicted to be metallic at the DFT level of theory, but 
semiconducting at the GW level. It thus provides an- 
other example besides bulk Ge where the DFT gap 
is not only quantitatively but also qualitatively wrong. 
Second, the difference between the DFT and QP gap can 
be understood in terms of the interaction between an ex- 
tra particle - which is missing in DFT - and the charge it 
induces on the molecular surface. Third, the Peierls dis- 
torted polyacene conformer is energetically favoured only 
for very long chains and posesses a QP spectrum which 
is markedly different from the aromatic form. This also 
underlines the usefulness of approximate GW schemes, 
since in a first principles context the Peierls transition 
found here might not be noticed due to the limited treat- 
able system size 62] . 



IV. NUMERICAL CONSIDERATIONS 

In the following, we shortly discuss the numerical ef- 
ficiency of our approximations. The method scales like 
N'^Nf, where N is the number of basis functions and 
Ni is the number of spherically averaged basisfunctions, 
that are used in the representation of two-point func- 
tions [Ni < N). This has to be compared to first princi- 
ples implementations, which usually scale like N*. An 
additional reduction of computation time is obtained, 
since we use a minimal basis of optimized atomic orbitals, 
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while in a first principles framework a larger number of 
primitive orbitals is required. Moreover, the scaling pref- 
actor is reduced in the DFTB scheme, because the nec- 
essary integrals are either precalculated and tabulated 
or approximated by simple functions. As an example, 
the first principles evaluation of the QP spectrum of an- 
thracene took 170 minutes on a Pentiimr Xeon 2.20 GHz 
processor (including 120 minutes for the DFT part of the 
calculation), compared to less than 1 second on a Pen- 
tium 4 with 2.40 GHz in our approach. The largest struc- 
ture we studied, the n = 30 polyacene with 186 atoms 
took 10 minutes. The limiting factor in the calculation 
of very large systems is therefore not the computation 
time but rather the memory requirement. We try to cir- 
cumvent this problem by computing memory intensive 
quantities like the overlap charges q^^ on-the-fly in a di- 
rect way. 

V. CONCLUDING REMARKS 

In this work we presented a method to perform quasi- 
particle calculations of molecular systems at a highly re- 
duced computational cost compared to first principles im- 



plementations. The scheme was applied to hydrocarbons, 
but it can be easily extended to other elements, since all 
required parameters are calculated from first principles. 
The various approximations of the method are intended 
to be as consistent as possible with the underlying DFTB 
approach to allow for a stable error cancellation. For 
benzene and anthracene the results are indeed compara- 
ble with higher level calculations with the exception of 
(T-orbitals. Here, ways to overcome the deficiencies were 
outlined. Nevertheless, we think that the scheme could 
be useful already at the present stage, since e.g. for op- 
tical spectra or in the electronic transport only a few 
states around the Fermi level are active and dominate 
the physical properties of a system. 
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